#ifndef LDPC3_H_INCLUDED
#define LDPC3_H_INCLUDED

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Gauss-Markov channel simulator for LDPC codes															      //
// "Genius" with memory
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//#include <windows.h>
//#include <stdio.h>
#include <time.h>
//#include <iostream>
#include <fstream>				//REMOVER!!!!!!!!!!!!!
//#include <math.h>
//#include <stdlib.h>
#include <vector>
//#include <C:\AMD\Support\LDPC\Headers\archives.h>

using namespace std;

int ldpc3(double mi, double pi1, double K, int iErrors, int iRounds, double Eb0, double EbN, double step, int iMax, int ene, int ca, double timeT, int threadNumb)
{
////////////////////////////////////DEFINE HERE YOUR SIMULATION PARAMETERS//////////////////////////////
    unsigned int n=ene,k=ca,itmax=iMax;     // n and k-> matrix sizes   itmax-> maximum number of iterations of Sum-Product Algorithm
	double EbNo=Eb0, EbNoend=EbN, stepsize=step;   //EbNos-> Start and end sizes of noise      stepsize-> size of increments
    double MI=mi, PI1=pi1;							//MI-> memory	PI1-> Probability Bad State
    double Kappa=K;								//A is the impulsive noise
////////////////////////////////////////////////////////////////////////////////////////////////////////

	double PI0;										//PI0-> Probability Good State
	double Q, q;
    unsigned int i,j,indicator=1,iterations=0,ecount1=0,rounds=1,errorcounter1=0;
    unsigned int gaux=0,g=0;
    unsigned int haux=0,h=0;

    double pi=3.14159265358979;
    double Lprod=1,Lsum=0,K1,K2;            //Sum-Product algorithm variables
    double arg1,arg2;
    double U1,U2,V1,V2,W,X,Y,Z;

//Time variables----------------------------------------------------------------

	clock_t sTime, fTime;

//Simulatiom parameters control variables---------------------------------------

	int nonZero = 0;
	char resume;

//H loader variables------------------------------------------------------------

	char ch;
    vector<vector<int> >H;                  //Input matrix
    vector<int> aux;

//Codeword, noise and variancy parameters---------------------------------------

	unsigned int state = 0;                 //State of channel
    vector<double> noise(n,0);
    vector<double> r(n,0);
    vector<double> f0(n,0),A1(n,0),sdVector(n,0);
	double EbNo_Pot,sdg,sdb,sd,BER;
	double R=(double)k/n;

//------------------------------------------------------------------------------

    vector<vector<double> >L1q(n-k,vector<double>(n,0)),L1r(n-k,vector<double>(n,0));
    vector<int> m1(n,0),S1(n-k,0),e1(n,0);   //Message and syndrome vectors

//Reloading last section--------------------------------------------------------

	if (threadNumb == 0)
	{
		ifstream outX("Results1.txt",ios::in | ios::binary);
		if(outX.fail()){}
		else
		{
			cout << "Thread " << (threadNumb + 1) << " - Resume last simulation? [y/n]" << endl;
			resume = getchar();
			if (resume == 'y')	outX >> iRounds >> iErrors >> Eb0 >> timeT >> mi >> pi1 >> K;
			else{}
			outX.close();
		}
	}

	if (threadNumb == 1)
	{
		ifstream outX2("Results2.txt",ios::in | ios::binary);
		if(outX2.fail()){}
		else
		{
			cout << "Thread " << (threadNumb + 1) << " - Resume last simulation? [y/n]" << endl;
			resume = getchar();
			if (resume == 'y')	outX2 >> iRounds >> iErrors >> Eb0 >> timeT >> mi >> pi1 >> K;
			else{}
			outX2.close();
		}
	}

	if (threadNumb == 2)
	{
		ifstream outX3("Results3.txt",ios::in | ios::binary);
		if(outX3.fail()){}
		else
		{
			cout << "Thread " << (threadNumb + 1) << " - Resume last simulation? [y/n]" << endl;
			resume = getchar();
			if (resume == 'y')	outX3 >> iRounds >> iErrors >> Eb0 >> timeT >> mi >> pi1 >> K;
			else{}
			outX3.close();
		}
	}

	if (threadNumb == 3)
	{
		ifstream outX4("Results4.txt",ios::in | ios::binary);
		if(outX4.fail()){}
		else
		{
			cout << "Thread " << (threadNumb + 1) << " - Resume last simulation? [y/n]" << endl;
			resume = getchar();
			if (resume == 'y')	outX4 >> iRounds >> iErrors >> Eb0 >> timeT >> mi >> pi1 >> K;
			else{}
			outX4.close();
		}
	}

//------------------------------------------------------------------------------
//Loading parity check matrix
//------------------------------------------------------------------------------
    ifstream in ("H.txt", ios::in |ios::binary);  //Parity check matrix

    if(!in)
    {
        cout<<"Cannot open file.\n";
        return 1;
    }

    while(in)
    {
        in.get(ch);


        if(in)  //Reading inicial pivotted matrix and loading vector
        {
            if(ch=='0')
            {
                aux.push_back(0);
            }
            if(ch=='1')
            {
                aux.push_back(1);
            }

            if(aux.size()==n)
            {
                H.push_back(aux);
                aux.clear();
            }
        }
    }
    in.close();
    aux.clear();

    cout<<"Parity check matrices loaded"<<endl;
//------------------------------------------------------------------------------
// Creating Index matrices for storing non zero positions of H matrices
//------------------------------------------------------------------------------
    for(i=0; i<n-k; i++)
    {
        for(j=0; j<n; j++)
        {
            if(H[i][j]==1)
            {
                gaux++;   //finding de 1s in the matrix
            }
            if(g<gaux)
            {
                g=gaux;   //storing the quantity of 1s in g
            }

        }
        gaux=0;
    }

    vector<vector<int> >Index(n-k, vector<int>(g,-1));
    gaux=0;

    for(i=0; i<n-k; i++)
    {
        for(j=0; j<n; j++)
        {
            if(H[i][j]==1)
            {
                Index[i][gaux]=j;       //putting the 1s in position in the Index matrix
                gaux++;                 //in the maximum column index (1s in lines)
            }
        }
        gaux=0;
    }

    for(j=0; j<n; j++)
    {
        for(i=0; i<n-k; i++)
        {
            if(H[i][j]==1)
            {
                haux++;   //finding the 1s in te matrix
            }
            if(h<haux)
            {
                h=haux;   //storing the quantity of 1s em h
            }

        }
        haux=0;
    }

    vector<vector<int> >Indexcol(h,vector<int>(n,-1) );
    haux=0;

    for(j=0; j<n; j++)
    {
        for(i=0; i<n-k; i++)
        {
            if(H[i][j]==1)
            {
                Indexcol[haux][j]=i;    //putting the 1s in position in the Index matrix
                haux++;                 //in the maximum line index (1s in columns)
            }
        }
        haux=0;
    }
    H.clear(); // This is to free the memory used for storing H.

    cout<<"Index matrices loaded"<<endl;
    cout<<"Parity check matrices cleared"<<endl;
    cout<<"Simulation started"<<endl;

//------------------------------------------------------------------------------
//Opening file for saving partial results
//------------------------------------------------------------------------------

    ofstream outP("Results1.txt",ios::in | ios::binary | ios::app );
	ofstream outP2("Results2.txt",ios::in | ios::binary | ios::app );
	ofstream outP3("Results3.txt",ios::in | ios::binary | ios::app );
	ofstream outP4("Results4.txt",ios::in | ios::binary | ios::app );

//------------------------------------------------------------------------------
	
	PI0 = 1 - PI1;									//PI for good state
	Q = PI1*(1 - MI);
	q = (1 - PI1)*(1 - MI);

	srand(time(NULL));
    EbNo=EbNo-stepsize;

	if (iRounds != 0)	nonZero = 1;

    while(EbNo<EbNoend)
    {
		sTime = clock();
        EbNo=EbNo+stepsize;

// Reinitializing constants---------------------------------------------------
        if (nonZero == 1)
		{
			rounds=iRounds;
			errorcounter1=iErrors;
			nonZero = 0;
		}
		else
		{
			timeT = 0;
			rounds=0;
			errorcounter1=0;
		}
        indicator=1;
        iterations=0;
        ecount1=0;
//-----------------------------------------------------------------------------

        EbNo_Pot = pow(10,(EbNo/10));		//Convert EbNo from dB
        sdg = 1.0/(sqrt(2*EbNo_Pot*R));		//Sigma (standard deviation) definition good
        sdb = sqrt(Kappa*sdg*sdg);				//Sigma bad

	Y = (double)rand()/RAND_MAX;				//Y=[0,1]

	if (Y < PI0)
	{  
		state = 0;								//Define the initial state Good
		sd  = sdg;
	}
	else
	{ 
		state = 1;								//Define the initial state Good
		sd = sdb;
	}

while(errorcounter1<100){ // Stops only when 100 bit errors are found.

//Implementing Gauss-Markov channel-----------------------------------------------------

for(int i = 0; i < n; i++)
{
	noise[i] = 0;
}
for(i=0;i<n;i++){

//Defining the state--------------------------------------------------------------------

Z = ((double)rand()/RAND_MAX)/10;				//Z=[0,1]

if(state == 0){								//for good state
	if(Z <= Q){
		state = 1;							//transition from good to bad
		sd = sdb;
	}
	if(Z > Q){
		state = 0;						//no transition
		sd = sdg;
	}
}
else if(state == 1) {										//for bad state
	if(Z <= q){
		state = 0;							//transition from bad to good
		sd = sdg;
	}
	if(Z > q){
		state = 1;						//no transition
		sd = sdb;
	}
}

sdVector[i] = sd;					//creating sd vector

//-----------------------------------------------------------------------------------------

   do
   {
	  U1= (double)rand()/RAND_MAX;           // U1=[0,1]
      U2= (double)rand()/RAND_MAX;           // U2=[0,1]
      V1=2*U1 -1;                            // V1=[-1,1]
      V2=2*U2 -1;                            // V2=[-1,1]
      W=V1*V1 + V2*V2;
   } while(W>=1);

   X=sqrt(-2 *log(W) / W) * V1;

   X = sd*X;     //Considering that the noise has zero mean

   noise[i] = X;

   //printf("%d\n", state);
}

int saida = 0;
for(i=0;i<n;i++){  // r is the received codeword
   r[i]=1+noise[i];
}

//Decoding----------------------------------------------------------------------

for(j=0;j<n;j++){ // Calculating 2*y/sd*sd
	f0[j]= 2*r[j]/sdVector[j]*sdVector[j];
}

            for(i=0; i<n-k; i++) //Initializing LLR matrix
            {
                for(j=0; j<g; j++)
                {
                    if(Index[i][j]!=-1)
                    {
                        L1q[i][Index[i][j]]=f0[Index[i][j]];
                    }
                    else
                    {
                        break;
                    }
                }
            }

//Sum Product Algorithm---------------------------------------------------------

            while(indicator==1 && iterations<itmax)
            {

//Calculating the Product-------------------------------------------------------

                for(i=0; i<(n-k); i++)                 //Update rule for check nodes
                {
                    for(j=0; j<g; j++)
                    {
                        if(Index[i][j]!=-1)
                        {
                            Lprod=Lprod*tanh(L1q[i][Index[i][j]]/2);     //Product
                        }
                        if(Index[i][j]==-1) break;
                    }

                    for(j=0; j<g; j++)
                    {
                        if(Index[i][j]!=-1)  //Note that atan(z)=1/2*[ln(1+z)-ln(1-z)]
                        {
                            K1= tanh(L1q[i][Index[i][j]]/2);
                            K2= (Lprod/K1);
                            arg1=1+K2;
                            if(arg1<=0)
                            {
                                arg1=1.0E-20;   // This conditional is to avoid underflow
                            }
                            arg2=1-K2;
                            if(arg2<=0)
                            {
                                arg2=1.0E-20;   // This conditional is to avoid underflow
                            }
                            L1r[i][Index[i][j]]=log(arg1)-log(arg2);     //Filling the matrix to next step
                        }
                        if(Index[i][j]==-1)break;
                    }
                    Lprod=1;
                }

//Calculating the Sum-----------------------------------------------------------

                for(j=0; j<n; j++)                //Update rule for variable nodes
                {
                    for(i=0; i<h; i++)
                    {
                        if(Indexcol[i][j]!=-1)
                        {
                            Lsum+=L1r[Indexcol[i][j]][j];       //Sum
                        }
                        if(Indexcol[i][j]==-1)break;
                    }

                    for(i=0; i<h; i++)
                    {
                        if(Indexcol[i][j]!=-1)
                        {
                            L1q[Indexcol[i][j]][j]=f0[j]+Lsum-L1r[Indexcol[i][j]][j];
                        }
                        if(Indexcol[i][j]==-1)break;
                    }
                    A1[j]=f0[j]+Lsum;         //Decision vector
                    Lsum=0;
                }

//Filling the message vector---------------------------------------------------

                for(j=0; j<n; j++)             //Termination Rule
                {
                    if(A1[j]>=0)
                    {
                        m1[j]=0;
                    }
                    else m1[j]=1;
                }

//------------------------------------------------------------------------------
                /*
                cout<<"The message vector 1 is: ";
                cout<<"\n";
                for(i=0;i<n;i++){
                        cout<<m1[i];
                }

                cout<<"\n";
                */




//Calculating the Syndrome------------------------------------------------------

                for(i=0; i<n-k; i++) //Initializing S[]
                {
                    S1[i]=0;
                }

                for(i=0; i<n-k; i++)   // Syndrome computation based on the Index matrix S= H*mt
                {
                    for(j=0; j<g; j++)
                    {
                        if(Index[i][j]!=-1)
                        {
                            S1[i]=(m1[Index[i][j]]+S1[i])%2;
                        }
                    }
                }


                for(j=0; j<n-k; j++)
                {
                    if(S1[j]!=0)
                    {
                        indicator=1;              //indicates if has Syndrome, if yes, repeats the Sum-Product Algorithm
                        break;
                    }
                    else indicator=0;         //if don't, step out to errors count
                }

//------------------------------------------------------------------------------

                iterations++;

            }

            indicator=1;
            iterations=0;

            ecount1=0;
            for(i=0; i<n; i++)
            {
                if(m1[i]!=0)
                {
                    ecount1++;   //message number of errors count
                }

            }

            if (ecount1>0)   // This indicates error in decoding
            {
                errorcounter1++;
            }

            rounds++;               //number of iterations


            if(rounds%1000==0)
            {
                cout<<"Round number: "<<rounds<<". Detected bit errors: "<<errorcounter1 << " Thread - " << (threadNumb + 1) <<endl;
				if(rounds%10000==0)
				{
					fTime = clock();
					
					if (threadNumb == 0)
					{
						//rounds, errors, EbNo, time
						outP.close();
						fopen("Results1.txt", "w"); 
						ofstream outP("Results1.txt",ios::in | ios::binary | ios::app ); //Time Log file of the simulation

						outP << rounds << " " << errorcounter1 << " " << EbNo << " " << (((double)(fTime - sTime) / CLOCKS_PER_SEC ) + timeT) << " " << MI << " " << PI1 << " " << Kappa;
						outP << (char)13 << (char)10;
					}

					if (threadNumb == 1)
					{
						//rounds, errors, EbNo, time
						outP2.close();
						fopen("Results2.txt", "w"); 
						ofstream outP2("Results2.txt",ios::in | ios::binary | ios::app ); //Time Log file of the simulation

						outP2 << rounds << " " << errorcounter1 << " " << EbNo << " " << (((double)(fTime - sTime) / CLOCKS_PER_SEC ) + timeT) << " " << MI << " " << PI1 << " " << Kappa;
						outP2 << (char)13 << (char)10;						
					}

					if (threadNumb == 2)
					{
						//rounds, errors, EbNo, time
						outP3.close();
						fopen("Results3.txt", "w"); 
						ofstream outP3("Results3.txt",ios::in | ios::binary | ios::app ); //Time Log file of the simulation

						outP3 << rounds << " " << errorcounter1 << " " << EbNo << " " << (((double)(fTime - sTime) / CLOCKS_PER_SEC ) + timeT) << " " << MI << " " << PI1 << " " << Kappa;
						outP3 << (char)13 << (char)10;						
					}

					if (threadNumb == 3)
					{
						//rounds, errors, EbNo, time
						outP4.close();
						fopen("Results4.txt", "w"); 
						ofstream outP4("Results4.txt",ios::in | ios::binary | ios::app ); //Time Log file of the simulation

						outP4 << rounds << " " << errorcounter1 << " " << EbNo << " " << (((double)(fTime - sTime) / CLOCKS_PER_SEC ) + timeT) << " " << MI << " " << PI1 << " " << Kappa;
						outP4 << (char)13 << (char)10;
					}
					
				}
            }

        }


        BER=(double)(errorcounter1)/(double)rounds;         //errors number/iterations number

		fTime = clock();

		ofstream ouTime("Times.txt",ios::in | ios::binary | ios::app ); //Time Log file of the simulation
		ofstream ouTime2("Times2.txt",ios::in | ios::binary | ios::app ); //Time Log file of the simulation
		ofstream ouTime3("Times3.txt",ios::in | ios::binary | ios::app ); //Time Log file of the simulation
		ofstream ouTime4("Times4.txt",ios::in | ios::binary | ios::app ); //Time Log file of the simulation

		if (threadNumb == 0)
		{
			ouTime << "EbNo=" << EbNo << " load time: " << (((double)(fTime - sTime) / CLOCKS_PER_SEC) + timeT) << " seconds.";
			ouTime << (char)13 << (char)10;	
		}

		if (threadNumb == 1)
		{
			ouTime2 << "EbNo=" << EbNo << " load time: " << (((double)(fTime - sTime) / CLOCKS_PER_SEC) + timeT) << " seconds.";
			ouTime2 << (char)13 << (char)10;
		}

		if (threadNumb == 2)
		{
			ouTime3 << "EbNo=" << EbNo << " load time: " << (((double)(fTime - sTime) / CLOCKS_PER_SEC) + timeT) << " seconds.";
			ouTime3 << (char)13 << (char)10;
		}

		if (threadNumb == 3)
		{
			ouTime4 << "EbNo=" << EbNo << " load time: " << (((double)(fTime - sTime) / CLOCKS_PER_SEC) + timeT) << " seconds.";
			ouTime4 << (char)13 << (char)10;
		}

		ouTime.close();
		ouTime2.close();
		ouTime3.close();
		ouTime4.close();

//------------------------------------------------------------------------------
//Opening file for saving final results
//------------------------------------------------------------------------------

        ofstream out("Thread1.txt",ios::in | ios::binary | ios::app ); // Log file of the simulation
        ofstream out2("Thread2.txt",ios::in | ios::binary | ios::app ); // Log file of the simulation
        ofstream out3("Thread3.txt",ios::in | ios::binary | ios::app ); // Log file of the simulation
        ofstream out4("Thread4.txt",ios::in | ios::binary | ios::app ); // Log file of the simulation

        //out << "Thread - " << (threadNumb + 1);
        //out<<(char)13<<(char)10;
        if (threadNumb == 0)
        {
            out<<"EbNo= "<<EbNo<<". Number of rounds for 100 block errors= "<<rounds;
            out<<(char)13<<(char)10;
            out<<"Total number of errors: "<<errorcounter1;
            out<<(char)13<<(char)10;
            out<<"BER: "<<BER;
            out<<(char)13<<(char)10;
            out<<(char)13<<(char)10;            //same as "endl" (multiplataform)
        }
        if (threadNumb == 1)
        {
            out2<<"EbNo= "<<EbNo<<". Number of rounds for 100 block errors= "<<rounds;
            out2<<(char)13<<(char)10;
            out2<<"Total number of errors: "<<errorcounter1;
            out2<<(char)13<<(char)10;
            out2<<"BER: "<<BER;
            out2<<(char)13<<(char)10;
            out2<<(char)13<<(char)10;            //same as "endl" (multiplataform)
        }
        if (threadNumb == 2)
        {
            out3<<"EbNo= "<<EbNo<<". Number of rounds for 100 block errors= "<<rounds;
            out3<<(char)13<<(char)10;
            out3<<"Total number of errors: "<<errorcounter1;
            out3<<(char)13<<(char)10;
            out3<<"BER: "<<BER;
            out3<<(char)13<<(char)10;
            out3<<(char)13<<(char)10;            //same as "endl" (multiplataform)
        }
        if (threadNumb == 3)
        {
            out4<<"EbNo= "<<EbNo<<". Number of rounds for 100 block errors= "<<rounds;
            out4<<(char)13<<(char)10;
            out4<<"Total number of errors: "<<errorcounter1;
            out4<<(char)13<<(char)10;
            out4<<"BER: "<<BER;
            out4<<(char)13<<(char)10;
            out4<<(char)13<<(char)10;            //same as "endl" (multiplataform)
        }


        out.close();
        out2.close();
        out3.close();
        out4.close();
		outP.close();
		outP2.close();
		outP3.close();
		outP4.close();

        //printing the same on screen
        cout << "Thread - " << (threadNumb + 1) << endl;
        cout<<"EbNo= "<<EbNo<<". Number of rounds for 100 block errors= "<<rounds<<endl;
        cout<<"Total number of errors: "<<errorcounter1<<endl;
        cout<<"BER "<<BER<<endl<<endl;
    }

    return 0;
}

#endif // LDPC3_H_INCLUDED
